#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "regularisation.h"
#include "spline.h"



int i;
void regular(double rhomx[][3], double cornergap, double dz , int iter_z){
    
    int id = 1;
    double uppergap = id+cornergap;
    double lowergap = id-cornergap;
    double gap_amount = uppergap-lowergap;
    int gapnode   = gap_amount/dz;

    double testvalue;
    int count_inter = 0;
    int leng  = (int) round((uppergap - lowergap)/dz-1);

    double spline[leng]; 
    int regcase;


    for(i = 0;i<iter_z;i++){
        testvalue = round(i*dz*100000)/100000;


        if ( (testvalue >= lowergap) && (testvalue <= uppergap)){



            if (  fabs(testvalue-lowergap)< pow(10,-7)){
                       splinefunc(spline,leng,lowergap,uppergap,dz);

                       if( (fabs(rhomx[i][0] - 1) <= pow(10,-7)) && (fabs(rhomx[i][1] - 0)<= pow(10,-7)) ){

                           // case 1
                           if ( (fabs(rhomx[i+gapnode][0] - 0)<=pow(10,-7)) && (fabs(rhomx[i+gapnode][1] - 1)<=pow(10,-7)) ){
                           regcase = 1;

                           }else{
                               //nothing case 1
                           }

                           // case 2
                           if ( (fabs(rhomx[i+gapnode][0] - 0)<=pow(10,-7)) && (fabs(rhomx[i+gapnode][2] - 1)<=pow(10,-7)) ){
                           regcase = 2;

                           }else{
                              //nothing case 2
                           }
                       }else{
                               //next
                       }
                       if( fabs(rhomx[i][0] - 0)<= pow(10,-7) && fabs(rhomx[i][1] - 1)<= pow(10,-7) ){
                           // case 3
                           if ( (fabs(rhomx[i+gapnode][0] - 1)<=pow(10,-7)) && (fabs(rhomx[i+gapnode][1] - 0)<=pow(10,-7)) ){
                               regcase = 3;

                           }else{
                               //nothing case 3
                           }
                           // case 4
                           if ( (fabs(rhomx[i+gapnode][2] - 1)<=pow(10,-7)) && (fabs(rhomx[i+gapnode][1] - 0)<=pow(10,-7)) ){
                               regcase = 4;

                           }else{
                               //noting case 4
                           }
                       }else{
                           //next
                       }
                       if( (rhomx[i][0] - 0)<= pow(10,-7) && (rhomx[i][1] - 0)<= pow(10,-7) ){
                           // case 5
                           if ( (fabs(rhomx[i+gapnode][0] - 1)<=pow(10,-7)) && (fabs(rhomx[i+gapnode][1] - 0)<=pow(10,-7)) ){
                               regcase = 5;
                           }else{
                               //nothing case 5
                           }
                           // case 6
                           if ( (fabs(rhomx[i+gapnode][0] - 1)<=pow(10,-7)) && (fabs(rhomx[i+gapnode][1] - 1)<=pow(10,-7)) ){
                               regcase = 6;
                           }else{
                               //nothing case 6
                           }
                       }else{
                           //next
                       }





        }else{



                if(testvalue > lowergap && testvalue < uppergap){

                switch (regcase){

                   case 0:
                   // no regularisation
                       break;
                   case 1:
                   rhomx[i][0] = spline[count_inter];
                   rhomx[i][1] = 1-spline[count_inter];
                       break;
                   case 2:
                   rhomx[i][0] = spline[count_inter];
                   rhomx[i][2] = 1-spline[count_inter];
                       break;
                   case 3:
                   rhomx[i][1] = spline[count_inter];
                   rhomx[i][0] = 1-spline[count_inter];
                       break;
                   case 4:
                   rhomx[i][1] = spline[count_inter];
                   rhomx[i][2] = 1-spline[count_inter];
                   break;
                   case 5:
                   rhomx[i][2] = spline[count_inter];
                   rhomx[i][0] = 1-spline[count_inter];
                   break;
                   case 6:
                   rhomx[i][2] = spline[count_inter];
                   rhomx[i][1] = 1-spline[count_inter];
                   break;
               default:
                 break;
                }//end switch

               count_inter++;

           }else{

           id++;
           uppergap = id+cornergap;
           lowergap = id - cornergap;
           regcase = 0;
           count_inter = 0;


                }
            }// test lowergap
   }else{

        }
    if (id == 3){
        break;
    }
  }
}//end reg

